for (package in c("tidyverse","fpp3", "GGally", "normtest")) {
if (!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.5 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.0.2 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## ── Attaching packages ──────────────────────────────────────────── fpp3 0.4.0 ──
## ✓ lubridate 1.8.0 ✓ feasts 0.2.2
## ✓ tsibble 1.1.1 ✓ fable 0.3.1
## ✓ tsibbledata 0.4.0
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x tsibble::intersect() masks base::intersect()
## x tsibble::interval() masks lubridate::interval()
## x dplyr::lag() masks stats::lag()
## x tsibble::setdiff() masks base::setdiff()
## x tsibble::union() masks base::union()
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
# Contraste para los coeficientes
my_t_test <- function (object, ...)
{
par <- rbind(t_stat=tidy(object)$statistic, p_value=tidy(object)$p.value)
colnames(par) <- tidy(object)$term
if (NCOL(par) > 0) {
cat("\nt-test:\n")
coef <- round(par, digits = 4)
print.default(coef, print.gap = 2)
}
}
# Gráfico de correlogramas de residuos
my_tsresiduals <- function (data, ...) {
if (!fabletools::is_mable(data)) {
abort("gg_tsresiduals() must be used with a mable containing only one model.")
}
data <- stats::residuals(data)
if (n_keys(data) > 1) {
abort("gg_tsresiduals() must be used with a mable containing only one model.")
}
gg_tsdisplay(data, !!sym(".resid"), plot_type = "partial",
...)
}
# Validación cruzada anidada
nested_cv <- function(df, h, last_train, string_formula){
nested_errors <- vector()
for (i in seq(last_train, last(df$date), h)){
train <- df %>% filter(date<=i)
test <- df %>% filter(date>i)
fitted_model <- train %>%
model(arima = ARIMA(as.formula(string_formula)))
h_forecast = min(dim(test)[1], h)
fc <- fitted_model %>%
forecast(h=h_forecast)
test_err <- fc %>%
accuracy(test) %>%
select(MAPE)
nested_errors <- c(nested_errors, test_err$MAPE)
NewList <- list("errors"=nested_errors, "mean"=mean(nested_errors))
}
return(NewList)}
# Test de autocorrelacines de los residuos
autocorrelation_test_plot <- function(aug, dof = 4, m = 12, h = 5, alpha = 0.05){
vec <- c()
num_lags = seq(1, h*m)
for (i in num_lags){
vec <- c(vec,aug %>% features(.resid, ljung_box, lag=i, dof=dof) %>% .$lb_pvalue)
}
autocorr_pvalues_resid <- tibble(
lag = num_lags,
p_value = vec,
incorelated = p_value >= alpha
)
plot <- autocorr_pvalues_resid %>%
drop_na() %>%
ggplot(aes(lag, p_value, color = incorelated)) +
geom_point() +
geom_hline(aes(yintercept = alpha), linetype="dashed", color = "indianred2")
newList <- list("values" = autocorr_pvalues_resid, "plot" = plot)
return(newList)
}
Antes de empezar a hacer el estudio hay que distinguir la parte de los datos que se utilizará para construir la fórmula (para modelizar) de aquella que se utilizará para validar los resultados y cuyos datos, no deben ser utilizados.
Convertimos la serie en un objeto tsibble, para trabajar de manera cómoda. Analizamos visualmente la serie y las marcas de tiempo para saber cómo dividir los datos.
TU <- read_csv('SerieTU_format.csv')
## Rows: 156 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (1): value
## date (1): fecha
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# arreglamos formato de fecha
TU$fecha = tsibble::yearmonth(TU$fecha)
TU = as_tsibble(TU, index = fecha)
TU %>%
autoplot(value) +
labs(title = "Transporte urbano") +
xlab("Fecha") + ylab("Miles de viajeros") +
geom_point(aes(y = value), size = 1.2, shape = 16, color = "black")
Dado que tenemos datos mensuales y visualmente se aprecia una estacionalidad anual, reservamos el último año para test y el resto para train.
TU_train <- TU %>% filter_index(. ~ "2012-12-01")
TU_test <- TU %>% filter_index("2013-01-01" ~ .)
Se comienza la fase de identificación donde se deciden las transformaciones a realizar y el propio ajuste de la serie. Aunque veremos que cuando se realiza la selección de los parámetros del modelo ARIMA, la iteración con la estimación y el contraste son constantes hasta converger a un modelo válido.
Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada.
TU_train %>%
gg_season(value, period = "year",labels = "left")
Estacionalidad anual: la forma de las series anuales son iguales.
Descomposición STL para analizar la existencia de tendencia y estacionalidad.
TU_train %>%
model(seats = feasts:::STL(value)) %>%
components() %>%
autoplot()
Se confirman visualmente la estacionalidad anual y una tendencia creciente y luego decreciente.
A continuación, estudiamos si la serie es estacionaria en media y varianza y actuamos en consecuencia.
Analicemos qué transformación Box-Cox puede aplicarse para estabilizar la varianza.
lambda <- TU_train %>%
features(value, features = guerrero) %>% pull(lambda_guerrero)
lambda
## [1] 1.876598
TU_train %>% autoplot(box_cox(value,lambda)) +
labs(y = "Box-Cox transformed TU")
Dado que la serie visualmente muestra una varianza más o menos constante, que el valor Box Cox es más próximo a 1 que a otros valores de transformaciones usuales (logaritmo, raíz cuadrada) y que no ha habido un cambio muy llamativo con la transformación, no realizaremos ningún cambio a priori. Si en el proceso de estimación no conseguimos un modelo adecuado, volveríamos a este punto para valorar de nuevo la transformación.
La serie muestra una tendencia creciente y luego decreciente, es decir, la media va variando a lo largo de la serie. Por tanto, es esperable que se necesite realizar al menos una diferencia para estabilizar la media. Formalmente, se comprueba mediante una prueba de raíces unitarias.
TU_train %>%
features(value, unitroot_kpss)
El p-valor es menor que 0.05, lo que indica que la hipótesis nula es rechazada. Es decir, hay evidencias estadísticas para rechazar la hipótesis nula de que la serie sea estacionaria. Por tanto, habría que realizar una diferencia regular.
TU_train %>%
features(difference(value, 1), unitroot_kpss)
Una vez realizada la diferencia regular, el p-valor del test indica que no hay evidencias estadísticas para rechazar la hipótesis de que la serie sea estacionaria en media.
Graficamos la serie diferenciada y ya observamos que el nivel de la serie no cambia:
TU_train %>% autoplot(difference(value, 1)) + labs(title = "Serie con una diferencia regular")
Debido a la estacionalidad de los datos, quizá deba realizarse una diferencia estacional. Para evaluarlo podemos recurrir a la función unitroot_nsdiffs(), esta evalúa el estadístico fuerza estacional \(F_S\) y sugiere una diferencia estacional si esta es mayor que 0.64.
TU_train %>%
mutate(turnover = difference(value,1)) %>%
features(turnover, unitroot_nsdiffs)
El método sugiere una diferencia estacional. En cualquier caso, para asegurar que la serie es estacionaria, en la fase de identificación del modelo debemos asegurarnos de que las raíces de la parte autoregresiva están fuera del círculo unitario (y, por tanto, 1/raíz dentro del círculo).
En los pasos previos se ha sugerido una diferencia regular y una estacional. Empezaremos teniendo en cuenta con la diferencia regular. Es decir, nos dirigimos a identificar los órdenes \(p, q, P, Q\) de un SARIMA (p,d,q)x(P,D,Q)\(_{12}\) donde ya sabemos que \(d=1\). Vemos que los residuos contienen mucha información.
fit <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(0,0,0, period = 12)))
fit %>% my_tsresiduals(lag_max = 36)
En la ACF se aprecia una estructura autoregresiva regular y estacional cuyos órdenes se confirman en la PACF. Comenzamos introduciendo \(P=1\).
Comencemos ajustando la parte estacional mediante un SARIMA (0,1,0)x(1,0,0)\(_{12}\).
fit2 <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,0) + PDQ(1,0,0, period=12)))
fit2 %>% my_tsresiduals(lag_max =36)
report(fit2)
## Series: value
## Model: ARIMA(0,1,0)(1,0,0)[12] w/ drift
##
## Coefficients:
## sar1 constant
## 0.8615 -95.9125
## s.e. 0.0367 942.5549
##
## sigma^2 estimated as 263289328: log likelihood=-1596.33
## AIC=3198.66 AICc=3198.83 BIC=3207.55
my_t_test(fit2)
##
## t-test:
## sar1 constant
## t_stat 23.502 -0.1018
## p_value 0.000 0.9191
Hemos comprobado que el parámetro de este modelo es significativo y no existen problemas de estacionariedad.
Introducimos ahora \(q=1\). Es decir, SARIMA (0,1,1)x(1,0,0)\(_{12}\).
fit3 <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(1,0,0, period=12)))
fit3 %>% my_tsresiduals(lag_max =36)
report(fit3)
## Series: value
## Model: ARIMA(0,1,1)(1,0,0)[12] w/ drift
##
## Coefficients:
## ma1 sar1 constant
## -0.8816 0.920 -21.7567
## s.e. 0.0428 0.025 64.3210
##
## sigma^2 estimated as 118838387: log likelihood=-1542.58
## AIC=3093.15 AICc=3093.44 BIC=3105
my_t_test(fit3)
##
## t-test:
## ma1 sar1 constant
## t_stat -20.578 36.8724 -0.3383
## p_value 0.000 0.0000 0.7357
Introducimos \(Q=1\). Es decir, SARIMA (0,1,1)x(1,0,1)\(_{12}\).
fit4 <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(1,0,1, period=12)))
fit4 %>% my_tsresiduals(lag_max =36)
report(fit4)
## Series: value
## Model: ARIMA(0,1,1)(1,0,1)[12] w/ drift
##
## Coefficients:
## ma1 sar1 sma1 constant
## -0.7892 0.9990 -0.8655 0.0141
## s.e. 0.0435 0.0017 0.1157 1.5368
##
## sigma^2 estimated as 75708833: log likelihood=-1520.22
## AIC=3050.44 AICc=3050.88 BIC=3065.26
my_t_test(fit4)
##
## t-test:
## ma1 sar1 sma1 constant
## t_stat -18.1272 574.3631 -7.4833 0.0091
## p_value 0.0000 0.0000 0.0000 0.9927
El coeficiente estimado para “sar1” se acerca peligrosamente a 1. Es necesario realizar una diferencia estacional.
Diferencia \(D=1\). Es decir, SARIMA (0,1,1)x(0,1,1)\(_{12}\).
fit5 <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(0,1,1) + PDQ(0,1,1, period=12)))
fit5 %>% my_tsresiduals(lag_max =36)
report(fit5)
## Series: value
## Model: ARIMA(0,1,1)(0,1,1)[12] w/ poly
##
## Coefficients:
## ma1 sma1 constant
## -0.8385 -0.9923 -115.6608
## s.e. 0.0448 0.1955 34.8462
##
## sigma^2 estimated as 64684904: log likelihood=-1377.47
## AIC=2762.94 AICc=2763.26 BIC=2774.44
my_t_test(fit5)
##
## t-test:
## ma1 sma1 constant
## t_stat -18.7347 -5.0763 -3.3192
## p_value 0.0000 0.0000 0.0012
Introducimos \(p=1\). Es decir, SARIMA (1,1,1)x(0,1,1)\(_{12}\).
fit6 <- TU_train %>%
model(arima = ARIMA(value ~ 1+pdq(1,1,1) + PDQ(0,1,1, period=12)))
fit6 %>% my_tsresiduals(lag_max =36)
report(fit6)
## Series: value
## Model: ARIMA(1,1,1)(0,1,1)[12] w/ poly
##
## Coefficients:
## ar1 ma1 sma1 constant
## -0.3627 -0.7216 -0.8861 -157.0759
## s.e. 0.0937 0.0657 0.1301 58.5297
##
## sigma^2 estimated as 64256797: log likelihood=-1371
## AIC=2751.99 AICc=2752.47 BIC=2766.37
my_t_test(fit6)
##
## t-test:
## ar1 ma1 sma1 constant
## t_stat -3.8717 -10.9771 -6.8126 -2.6837
## p_value 0.0002 0.0000 0.0000 0.0082
Introducimos \(p=2\) y quitamos constante. Es decir, SARIMA (2,1,1)x(0,1,1)\(_{12}\).
fit7 <- TU_train %>%
model(arima = ARIMA(value ~ pdq(2,1,1) + PDQ(0,1,1, period=12)))
fit7 %>% my_tsresiduals(lag_max =36)
report(fit7)
## Series: value
## Model: ARIMA(2,1,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sma1
## -0.7184 -0.3967 -0.3891 -0.7876
## s.e. 0.1308 0.1117 0.1412 0.0879
##
## sigma^2 estimated as 64713127: log likelihood=-1368.38
## AIC=2746.76 AICc=2747.24 BIC=2761.13
my_t_test(fit7)
##
## t-test:
## ar1 ar2 ma1 sma1
## t_stat -5.4916 -3.5525 -2.7562 -8.9572
## p_value 0.0000 0.0005 0.0067 0.0000
gg_arma(fit7)
Para completar la fase de contraste evaluamos los residuos. Veamos primero el histograma.
aug <-fit7 %>% augment()
# Histogram
aug %>%
ggplot(aes(x = .resid)) +
geom_histogram(bins = 50) +
ggtitle("Histogram of residuals")
Contraste t-student para contrastar si la media es 0.
# Student's t-Test for mean=0
t.test(aug$.resid)
##
## One Sample t-test
##
## data: aug$.resid
## t = -0.69101, df = 143, p-value = 0.4907
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -1682.5268 810.8847
## sample estimates:
## mean of x
## -435.8211
Como p-value > 0.05 no existen evidencias significativas para rechazar la hipótesis nula de que la muestra tenga media 0.
A continuación comprobamos que los residuos están incorrelados.
# Ljung-Box autocorrelation
aug %>% features(.resid, ljung_box, lag=12, dof=4)
Como el p-valor es >0.05, podemos concluir que los residuos están incorrelados.
Para contrastar la heterocedasticidad se puede utilizar una regresión media-dispersión. Calculando los grupos de manera anual, dado que esa es la estacionalidad de la serie.
log_log <- aug %>% as_tibble() %>%
group_by(year(fecha)) %>%
summarize(mean_resid = log(mean(.resid+1)), std_resid = log(sd(.resid+1)))
summary(lm(std_resid~mean_resid, log_log))
##
## Call:
## lm(formula = std_resid ~ mean_resid, data = log_log)
##
## Residuals:
## 1 2 3 5
## 0.05104 0.10698 -0.51048 0.35246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.319 0.718 1.837 0.2076
## mean_resid 0.986 0.106 9.304 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4466 on 2 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.9774, Adjusted R-squared: 0.9661
## F-statistic: 86.56 on 1 and 2 DF, p-value: 0.01136
El p-valor de log(media) es menor que 0.05, se rechaza la hipótesis de que sean homocedásticos.
Para contrastar la normalidad realizamos el test Jarque-Bera para evaluar la normalidad.
# Jarque Bera test
jb.norm.test(na.omit(aug$.resid))
##
## Jarque-Bera test for normality
##
## data: na.omit(aug$.resid)
## JB = 32.549, p-value = 0.001
Como el p-valor es menor que el nivel de significancia 0.05, se rechaza la hipótesis nula de normalidad de los residuos.
Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones en el conjunto de test, calculamos sus errores y los comparamos con los residuos.
# Residual accuracy
resids <- fit7 %>%
accuracy() %>%
select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>%
mutate(Evaluation='Training')
# Forecasting
fc <- fit7 %>%
forecast(h=12)
test_err <- fc %>%
accuracy(TU) %>%
select(-c(.model, .type, ME, MPE, ACF1, RMSSE)) %>%
mutate(Evaluation='Test')
# Show errors together
bind_rows(test_err, resids) %>% select(Evaluation, everything())
Además de evaluar los errores con estas medidas globales, podemos evaluar los errores cometidos año a año visualmente. Así, podemos detectar outliers o efectos de calendario no detectados. Pudiendo así corregirlos con modelos más complejos mediante variables exógenas.
aug %>% ggplot() +
geom_line(aes(x = fecha, y = .fitted), color="navy") +
geom_line(aes(x = fecha, y = value), color="gray24") +
# geom_line(data=demandaGas_forecast, aes(x = index, y = value), color="red4") +
ggtitle("SARIMA train fitted values") +
xlab('Dates') +
ylab('Demand') + facet_wrap(vars(year(fecha)), scales = 'free')
También podemos dibujar las predicciones con los intervalos de confianza (aunque en este caso no sean representativos, dada la no normalidad de los residuos). Evaluamos así también la capacidad de generalización del modelo, comparando las predicciones con el conjunto de test.
h <- dim(TU_test)[1]
demanda_plot <- TU %>% filter(fecha>last(TU_train$fecha)-14 & fecha< last(TU_train$fecha) + 12 )
fit7 %>%
forecast(h=12) %>%
autoplot(demanda_plot)
Como sabemos que los modelos SARIMA no predicen bien más alla del orden más alto disponible, en este caso 12, los errores en el conjunto de test más allá de 12 datos se dipararían.